cleanup before you start
rm( list=ls() )
source('../lib/utilities.R')
library(moments)
save.plots = TRUE
setup
D <- read.table('../data/gdp-us-grate.csv')
dates <- as.Date(as.character(D[,1]),'%Y-%m-%d')
H <- 6
N <- nrow(D)-6
y <- D[1:N,2]
y.out <-D[(N+1):(N+H),2]
# LBQ Test
Box.test( y, lag=22 , type="Ljung-Box" )
##
## Box-Ljung test
##
## data: y
## X-squared = 84.824, df = 22, p-value = 2.563e-09
# ACF & PACF
par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( y , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( y , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)
Estimate ARMA Models
ar1 = arima(y,order=c(1,0,0))
ma1 = arima(y,order=c(0,0,1))
ma2 = arima(y,order=c(0,0,2))
arma11 = arima(y,order=c(1,0,1))
# GoF
ar1_aic <- (-2*ar1$loglik+2*3)/N
ma1_aic <- (-2*ma1$loglik+2*3)/N
ma2_aic <- (-2*ma2$loglik+2*4)/N
arma11_aic <- (-2*arma11$loglik+2*4)/N
ar1_bic <- (-2*ar1$loglik+log(N)*3)/N
ma1_bic <- (-2*ma1$loglik+log(N)*3)/N
ma2_bic <- (-2*ma2$loglik+log(N)*4)/N
arma11_bic <- (-2*arma11$loglik+log(N)*4)/N
round( rbind( c(ar1$loglik,ma1$loglik.ma2$loglik,arma11$loglik),
c(ar1_aic,ma1_aic,ma2_aic,arma11_aic) ,
c(ar1_bic,ma1_bic,ma2_bic,arma11_bic) ) , 3 )
## [,1] [,2] [,3] [,4]
## [1,] -767.611 -766.537 -767.611 -766.537
## [2,] 5.351 5.389 5.341 5.351
## [3,] 5.390 5.427 5.391 5.402
FITTED VALUES
ar1_mu <- y-ar1$residuals
ma1_mu <- y-ma1$residuals
ma2_mu <- y-ma2$residuals
arma11_mu <- y-arma11$residuals
ar1_res <- as.numeric(ar1$residuals)
ma1_res <- as.numeric(ma1$residuals)
ma2_res <- as.numeric(ma2$residuals)
arma11_res <- as.numeric(arma11$residuals)
Plot a bunch of stuff..
# AR1
par( mar=c(2,2,1,1) , xaxs="i" , mfrow=c(2,1) )
myplot( dates[1:N] , ar1_mu , t='l' , lwd=2 , col='blue3' , ylim=c(0,10) )
grid( lwd=1 , col="darkgrey" )
myplot( dates[1:N] , ar1_res/sd(ar1_res) , col='purple' , t='p' , ylim=c(-4,4) )
abline( h=0 , lwd=2 )
grid( lwd=1 , col="darkgrey" )
par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( ar1_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( ar1_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)
# MA1
par( mar=c(2,2,1,1) , xaxs="i" , mfrow=c(2,1) )
myplot( dates[1:N] , ma1_mu , col='blue3' , ylim=c(0,10) )
grid( lwd=1 , col="darkgrey" )
myplot( dates[1:N] , ma1_res/sd(ma1_res) , col='purple' , t='p' , ylim=c(-4,4) )
abline( h=0 , lwd=2 )
grid( lwd=1 , col="darkgrey" )
par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( ma1_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( ma1_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)
# MA2
par( mar=c(2,2,1,1) , xaxs="i" , mfrow=c(2,1) )
myplot( dates[1:N] , ma2_mu , t='l' , lwd=2 , col='blue3', ylim=c(0,10) )
grid( lwd=1 , col="darkgrey" )
myplot( dates[1:N] , ma2_res/sd(ma2_res) , col='purple' , t='p' , ylim=c(-4,4) )
abline( h=0 , lwd=2 )
grid( lwd=1 , col="darkgrey" )
par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( ma2_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( ma2_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)
# ARMA11
par( mar=c(2,2,1,1) , xaxs="i" , mfrow=c(2,1) )
myplot( dates[1:N] , arma11_mu , lwd=2 , col='blue3' , ylim=c(0,10))
grid( lwd=1 , col="darkgrey" )
myplot( dates[1:N] , arma11_res/sd(arma11_res) , col='purple' , t='p' , ylim=c(-4,4) )
abline( h=0 , lwd=2 )
grid( lwd=1 , col="darkgrey" )
par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( arma11_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( arma11_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)
Residuals
par( mar=c(2,2,3,2) , mfrow=c(2,2) )
kernel <- density(ar1_res/sqrt(ar1$sigma2))
plot( kernel , main='AR1' )
polygon( kernel , col="tomato" , border='darkred')
abline(h=0,lwd=2)
lines( seq(-10,20,0.1) , dnorm( seq(-10,20,0.1) ) , col='darkblue' ,lwd=2 )
kernel <- density(ma1_res/sqrt(ma1$sigma2))
plot( kernel , main='MA1' )
polygon( kernel , col="tomato" , border='darkred')
abline(h=0,lwd=2)
lines( seq(-10,20,0.1) , dnorm( seq(-10,20,0.1) ) , col='darkblue' ,lwd=2 )
kernel <- density(ma2_res/sqrt(ma2$sigma2))
plot( kernel , main='MA2' )
polygon( kernel , col="tomato" , border='darkred')
abline(h=0,lwd=2)
lines( seq(-10,20,0.1) , dnorm( seq(-10,20,0.1) ) , col='darkblue' ,lwd=2 )
kernel <- density(arma11_res/sqrt(arma11$sigma2))
plot( kernel , main='ARMA11' )
polygon( kernel , col="tomato" , border='darkred')
abline(h=0,lwd=2)
lines( seq(-10,20,0.1) , dnorm( seq(-10,20,0.1) ) , col='darkblue' ,lwd=2 )
par( mar=c(2,2,3,2) , mfrow=c(2,2) )
qqnorm(ar1_res,col='tomato',main='AR1')
qqline(ar1_res,lwd=2,lty=3)
qqnorm(ma1_res,col='tomato',main='MA1')
qqline(ma1_res,lwd=2,lty=3)
qqnorm(ma2_res,col='tomato',main='MA2')
qqline(ma2_res,lwd=2,lty=3)
qqnorm(arma11_res,col='tomato',main='ARMA11')
qqline(arma11_res,lwd=2,lty=3)
FORECAST EXERCISE 1
ar1_pred <- predict( ar1 , n.ahead=H )
ma1_pred <- predict( ma1 , n.ahead=H )
ma2_pred <- predict( ma2 , n.ahead=H )
arma11_pred <- predict( arma11 , n.ahead=H )
ar1_mse = sqrt( mean( (y.out - as.numeric(ar1_pred$pred) )**2 ) )
ma1_mse = sqrt( mean( (y.out - as.numeric(ma1_pred$pred) )**2 ) )
ma2_mse = sqrt( mean( (y.out - as.numeric(ma2_pred$pred) )**2 ) )
arma11_mse = sqrt( mean( (y.out - as.numeric(arma11_pred$pred) )**2 ) )
myplot( dates[(N-10):(N+H)] , c(y[(N-10):N], y.out) , t='b', main=sprintf('AR(1) RMSE %3.3f',ar1_mse) , ylim=c(0,7) , col='darkorange' )
abline( v=dates[N] , lwd=2 )
abline( h=ar1$coef['intercept'] , lwd=2 )
lines( dates[(N-10):N] , ar1_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( dates[(N+1):(N+H)] , as.numeric(ar1_pred$pred) , t='b' , lwd=2 , col='blue3' )
myplot( dates[(N-10):(N+H)] , c(y[(N-10):N], y.out) , t='b' , main=sprintf('MA(1) RMSE %3.3f',ma1_mse) , ylim=c(0,7) , col='darkorange' )
abline( v=dates[N] , lwd=2 )
abline( h=ma1$coef['intercept'] , lwd=2 )
lines( dates[(N-10):N] , ma1_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( dates[(N+1):(N+H)] , as.numeric(ma1_pred$pred) , t='b' , lwd=2 , col='blue3' )
myplot( dates[(N-10):(N+H)] , c(y[(N-10):N], y.out) , t='b', main=sprintf('MA(2) RMSE %3.3f',ma2_mse) , ylim=c(0,7) , col='darkorange' )
abline( v=dates[N] , lwd=2 )
abline( h=ma2$coef['intercept'] , lwd=2 )
lines( dates[(N-10):N] , ma2_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( dates[(N+1):(N+H)] , as.numeric(ma2_pred$pred) , t='b' , lwd=2 , col='blue3' )
myplot( dates[(N-10):(N+H)] , c(y[(N-10):N], y.out) , t='b' , main=sprintf('ARMA(1,1) RMSE %3.3f',arma11_mse) , ylim=c(0,7) , col='darkorange' )
abline( v=dates[N] , lwd=2 )
abline( h=arma11$coef['intercept'] , lwd=2 )
lines( dates[(N-10):N] , arma11_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( dates[(N+1):(N+H)] , as.numeric(arma11_pred$pred) , t='b' , lwd=2 , col='blue3' )
round(c( ar1_mse , ma1_mse , ma2_mse , arma11_mse )*100 , 3 )
## [1] 1997.261 1997.538 1997.163 1996.688
FORECAST EXERCISE 2
ar1_pred <- rep(0,H)
ma1_pred <- rep(0,H)
ma2_pred <- rep(0,H)
arma11_pred <- rep(0,H)
for( m in 0:(H-1) ){
y <- D[1:(N+m),2]
ar1 = arima(y,order=c(1,0,0))
ma1 = arima(y,order=c(0,0,1))
ma2 = arima(y,order=c(0,0,2))
arma11 = arima(y,order=c(1,0,1))
ar1_pred[1+m] <- predict( ar1 , n.ahead=1 )$pred
ma1_pred[1+m] <- predict( ma1 , n.ahead=1 )$pred
ma2_pred[1+m] <- predict( ma2 , n.ahead=1 )$pred
arma11_pred[1+m] <- predict( arma11 , n.ahead=1 )$pred
}
ar1_r2 = 1-mean( (y.out - ar1_pred )**2 )/mean( (mean(y)-y.out)**2 )
ma1_r2 = 1-mean( (y.out - ma1_pred )**2 )/mean( (mean(y)-y.out)**2 )
ma2_r2 = 1-mean( (y.out - ma2_pred )**2 )/mean( (mean(y)-y.out)**2 )
arma11_r2 = 1-mean( (y.out - arma11_pred )**2 )/mean( (mean(y)-y.out)**2 )
myplot( dates[(N-10):(N+H)] , c(y[(N-10):N], y.out) , t='b', main=sprintf('AR(1) R2 %3.3f',ar1_r2) , ylim=c(0,7) , col='darkorange' )
abline( v=dates[N] , lwd=2 )
abline( h=ar1$coef['intercept'] , lwd=2 )
lines( dates[(N-10):N] , ar1_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( dates[(N+1):(N+H)] , ar1_pred , t='b' , lwd=2 , col='blue3' )
myplot( dates[(N-10):(N+H)] , c(y[(N-10):N], y.out) , t='b' , main=sprintf('MA(1) R2 %3.3f',ma1_r2) , ylim=c(0,7) , col='darkorange' )
abline( v=dates[N] , lwd=2 )
abline( h=ma1$coef['intercept'] , lwd=2 )
lines( dates[(N-10):N] , ma1_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( dates[(N+1):(N+H)] , ma1_pred , t='b' , lwd=2 , col='blue3' )
myplot( dates[(N-10):(N+H)] , c(y[(N-10):N], y.out) , t='b', main=sprintf('MA(2) R2 %3.3f',ma2_r2) , ylim=c(0,7) , col='darkorange' )
abline( v=dates[N] , lwd=2 )
abline( h=ma2$coef['intercept'] , lwd=2 )
lines( dates[(N-10):N] , ma2_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( dates[(N+1):(N+H)] , ma2_pred , t='b' , lwd=2 , col='blue3' )
myplot( dates[(N-10):(N+H)] , c(y[(N-10):N], y.out) , t='b' , main=sprintf('ARMA(1,1) R2 %3.3f',arma11_r2) , ylim=c(0,7) , col='darkorange' )
abline( v=dates[N] , lwd=2 )
abline( h=arma11$coef['intercept'] , lwd=2 )
lines( dates[(N-10):N] , arma11_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( dates[(N+1):(N+H)] , arma11_pred , t='b' , lwd=2 , col='blue3' )
round(c( ar1_r2 , ma1_r2 , ma2_r2 , arma11_r2 )*100 , 3 )
## [1] -42.544 -25.992 -39.463 -41.733